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In this paper, spatially high-order Residual-Distribution (RD) schemes using the first- 
order hyperbolic system method are proposed for general time-dependent advection-dif- 
fusion problems. The corresponding second-order time-dependent hyperbolic advection- 
diffusion scheme was first introduced in [NASA/TM-2014-218175, 2014], where rapid con- 
vergences over each physical time step, with typically less than five Newton iterations, were 
shown. In that method, the time-dependent hyperbolic advection-diffusion system (linear 
and nonlinear) was discretized by the second-order upwind RD scheme in a unified manner, 
and the system of implicit-residual-equations was solved efficiently by Newton’s method 
over every physical time step. In this paper, two techniques for the source term discretiza- 
tion are proposed; 1) reformulation of the source terms with their divergence forms, and 
2) correction to the trapezoidal rule for the source term discretization. Third-, fourth, and 
sixth-order RD schemes are then proposed with the above techniques that, relative to the 
second-order RD scheme, only cost the evaluation of either the first derivative or both the 
first and the second derivatives of the source terms. A special fourth-order RD scheme is 
also proposed that is even less computationally expensive than the third-order RD schemes. 

The second-order Jacobian formulation was used for all the proposed high-order schemes. 

The numerical results are then presented for both steady and time-dependent linear and 
nonlinear advection-diffusion problems. It is shown that these newly developed high-order 
RD schemes are remarkably efficient and capable of producing the solutions and the gra- 
dients to the same order of accuracy of the proposed RD schemes with rapid convergence 
over each physical time step, typically less than ten Newton iterations. 

I. Introduction 

In this paper, we construct very efficient high-order schemes for general time-dependent advection- 
diffusion problems, based on the first-order hyperbolic system method. ’ In this method, the diffusion term 
is reformulated as a hyperbolic system, leading to the unification of advection and diffusion as a single hyper- 
bolic system. The drastic change in the type of equations, parabolic to hyperbolic, brings several dramatic 
improvements in the construction of numerical schemes: hyperbolic schemes for diffusion, the same order of 
accuracy for the solution and the gradients, orders-of-magnitude convergence acceleration, etc., which have 
been demonstrated for steady diffusion and viscous problems in Refs. “ and unsteady advection-diffusion 
problems in Ref. * 1 It is based on the reformulation of the governing equations, and therefore applicable to 
any discretization method. In this work, we consider a Residual-Distribution (RD) method, which has been 
well developed for hyperbolic systems and has a superior feature of achieving second-order accuracy in a 
compact stencil. 

In the previous work, we extended the hyperbolic method, for the first time, to time-accurate com- 
putations by an implicit time-integration method based on the second-order backward difference formula. 
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The resulting scheme was applied to various time-dependent problems, demonstrating second-order accuracy 
for the solution and the gradient achieved at all interior and boundary nodes in uniform and nonuniform 
grids at every physical time step, and rapid convergence for solving implicit-residual equations by Newton’s 
method (i.e., less than 5 iterations per physical time step), which is possible by the compactness of the RD 
schemes. As a consequence of the first-order re-formulation of the equation, the number of linear relaxations 
performed at every Newton iteration was shown to increase only linearly with the grid size, not quadratically 
as typical for diffusion problems. The efficiency of the developed second-order schemes was demonstrated 
for linear and nonlinear advection-diffusion problems on highly refined grids, up to 30000 nodes. 

In this paper, we propose a very simple extension of the second-order schemes to higher-order. We 
show that high-order spatial accuracy can be achieved simply by modifying the source term discretization. 
There are two approaches to the source term discretization: 1) reformulation of the source terms with their 
divergence forms, and 2) correction to the trapezoidal rule for the source term discretization. The former 
technique is based on the divergence formulation of source terms proposed in Ref.: write the source term 

in the divergence form and discretize it in the same way as the flux divergence term. The latter is based 
on a high-order correction to the trapezoidal rule, and thus called here the generalized trapezoidal rule. 
In either case, high-order accuracy is achieved by making low-order truncation error terms proportional to 
the residual, which thus vanish in the steady state and yield high-order accuracy. We solve the resulting 
implicit-residual equations by an implicit solver based on the second-order Jacobian matrix developed in the 
previous work. As we will show, the implicit solver is as powerful as Newton’s method; e.g., eight orders 
of magnitude reduction can be achieved in 10 iterations. To enable time-accurate computations, we employ 
high-order versions of the backward difference formulas (BDF), which are incorporated as source terms, and 
solve the implicit-residual equations by the implicit solver over each physical time step. In this manner, the 
steady state is made equivalent to the next physical time with all the benefits of the hyperbolic method 
retained. 

The high-order RD schemes developed in this work are significantly different from other high-order RD 
schemes in that our schemes are based on the first-order hyperbolic system formulation of the advection- 
diffusion equation. In this approach, the loss of high-order accuracy in the intermediate Reynolds number, 
as discussed in Refs.,' 1- cannot occur because the advective and diffusive terms are fully integrated into a 
single hyperbolic system. If the original advection-diffusion equation is discretized, a high-order RD scheme 
needs to be developed for the diffusion term (i.e., second derivative) and then carefully combined with an 
advection scheme, e.g. by using a blending parameter as described in Ref., to avoid the loss of accuracy. 
Furthermore, while high-order RD schemes based on high-order elements require extra degrees of freedom for 
each variable, our schemes are based on linear elements for any order of accuracy but require extra gradient 
variable to be added to the solution vector. Note that the number of extra variables in the high-order elements 
increase for higher-order accuracy, but the number of extra variables required in our approach is fixed and 
independent of the order of accuracy. Our approach is somewhat similar to those in Refs., 1J- but again 
is significantly different by the use of first-order hyperbolic system formulation of the advection-diffusion 
equation and by the source term discretization techniques. It is emphasized that our schemes require only 
the first and second derivatives of the source term, or in some cases the first derivatives only; they do not 
require the gradient computation for the solution variables. 

The third-order schemes developed in this paper are similar to the third-order finite-volume scheme of 
Katz and Sankaran ■ in that the second-order truncation error is eliminated by making it proportional to 
the residual and the upgrade is achieved by second-order accurate gradients. However, as we demonstrate 
in this paper, the proposed high-order RD schemes have several superior features: 1) implicit solver can be 
constructed by the Jacobian derived from a compact second-order RD scheme, 2) gradient computations 
are required for the source terms only, and not for the solution, 3) stiffness due to the second derivative of 
the diffusion term is completely eliminated, 4) higher-order schemes can be constructed beyond third-order, 
5) the same order of accuracy is achieved for the gradients, as well. In particular, the fourth-order scheme 
constructed in Section 5 is remarkably more efficient because it does not require second derivatives of the 
source term, which are required in the schemes described in Refs. ■ 

In this paper, we focus on one-dimensional linear and nonlinear advection-diffusion problems. It cer- 
tainly serves as a basis for the development of high-order multi-dimensional RD schemes for more complex 
equations. Yet, more importantly, the one-dimensional high-order schemes developed in this paper could po- 
tentially bring significant improvements to practical problems such as material thermal response calculations 
of thermal protection systems of atmospheric entry vehicles, 1 ' -1 and the experimental aeroheating data re- 
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duction,- ’ 1 which are based on one-dimensional analyses and routinely used in industries (e.g. NASA). The 
extension to higher dimensions is beyond the scope of the paper; it will be addressed in a subsequent paper. 

The paper is organized as follows. In the next section, the time-dependent hyperbolic advection-diffusion 
system is described. In Section 3, a compact second-order residual-distribution scheme, a steady solver, and 
the second-order discretization are discussed. In Section 4, the third- and fourth-order RD schemes with 
source term reformulation are proposed. In Section 5, the third-, fourth, and sixth-order RD schemes with 
source term discretization are developed and proposed. Numerical results are then presented in Section 6. 
Finally, Section 7 concludes the study with remarks. 


II. Time-dependent Hyperbolic Advection-Diffusion System 


We start with a linear advection-diffusion equation to simplify the discussion. We will extend the discus- 
sion later to a more general nonlinear advection-diffusion equation. 

Consider the one-dimensional (1-D) time-dependent advection-diffusion equation: 

d t u + ad x u = isd xx u + S(x), (1) 

where a and is are both taken to be positive constant, and S is the source term. We will follow the procedure 
we described in Ref. and rewrite the above equation as a first-order hyperbolic advection-diffusion system: 

a 

d T u = — a d x u + is d x p — — u + S(x), (2) 

d T p = ( d x u-p)/T r , (3) 


where the relaxation time, T r > 0, is arbitrary and defined as described in Ref., and S includes any existing 
source terms from the advection-diffusion problem, S, as well as any additional terms that arise from the 
implicit time-stepping scheme, At is the physical time steps, and r is the pseudo time step. Note that the 
dtp is taken as pseudo time derivative, d T p. 

The variable a depends on the order of the Backward-Differencing-Formula (BDF); 1 for the lst-order 
(BDF1), 3/2 for the second-order (BDF2), 11/6 for the third-order (BDF3), 25/12 for the fourth-order, 
and 147/60 for the sixth-order time discretizations. The remaining terms in the BDF are stored in the 
source term function S{x). Towards the pseudo steady state state, the variable p approaches the solution 
gradient and hence the above equation becomes identical to the original advection-diffusion equation with 
the time derivative discretized by the BDF, i.e., a consistent discretization of the original equation. The 
system reduces, of course, to the original steady equation also in the physical steady state when d t u = 0. 
Second-order accurate unsteady computations have been demonstrated based on the above formulation in 
Ref. 6 

In the vector form, our time-dependent first-order advection-diffusion system can be written as 


where 


<9U 

Ifr 


+ A 


<9U 

dx 


= S, 


( 4 ) 


u = 

U 

, A = 

a —v 

, S = 

— au/At + S(x) 


V 


— 1/T r 0 


-p/T r 


For any positive T r , the Jacobian has the following two real eigenvales: 


Ai = 


1-W1 + 


4v 

a?T r 


Aa — 


1 + 1/1 


4 v 
a^fr 


( 5 ) 

( 6 ) 


with linearly independent eigenvectors, and therefore the above system is hyperbolic in the pseudo time. 
The hyperbolicity serves mainly as a guide for discretization: various discretization techniques are available 
for hyperbolic systems, e.g., upwinding. In addition to the convenience in discretization, the major benefits 
are: 1) the hyperbolic discretization results in a strong coupling between the two variables that results in the 
equal order of accuracy for both u and p = d x u on arbitrary grids, and 2) 0(l/h) acceleration is achieved 
in iterative convergence for the linearized residual equation in implicit solvers due to the 0(l/h) condition 
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number of the Jacobian, which is 0(h) smaller than that of the Jacobian derived from a conventional 
diffusion scheme. For explicit schemes, the 0(1 /h) acceleration is achieved by the 0(h) time step typical in 
methods for hyperbolic systems, which is 0(1 /h) times larger than a typical time step of 0(h 2 ) for diffusion. 
The 0(l/h) acceleration in convergence has been demonstrated over traditional methods for the diffusion 
equation, ! ’ for the advection-diffusion equation, 2 ’ for the compressible Navier-Stokes equations, and most 
recently for time-dependent linear and nonlinear advection-diffusion problems using RD technique. 

In the next sections, we first briefly describe the second-order time-dependent discretization process using 
the RD scheme. Further details are given in Ref. 1 ' We then extend the order of accuracy of the scheme to 
third-order and fourth-order RD hyperbolic advection-diffusion schemes with reformulation of the hyperbolic 
system, and finally to higher order (fourth- and sixth-order) RD hyperbolic advection-diffusion schemes with 
introduction of new source term discretization technique. 


III. Second-Order RD Hyperbolic Advection-Diffusion 

The RD method requires 1) evaluation of the cell (or element) residuals, and 2) the distribution of the 
residuals to the nodes bounding the cell. Consider a one-dimensional domain discretized with N nodes 
that are distributed arbitrary over the domain of interest with the solution, u, and the solution gradient, 
p = d x ii, data stored at the nodes denoted by a *=1,2, 3, ..., N (Fig.l). We define the cell-residual <& c by 

h L h R 

o o o 

i-1 i i+1 


Figure 1. Schematic of grid spacing for a 1-D grid. 


integrating the spatial part of the Eq. (4) over the cell, C , defined by the nodes * and * + 1: 





(-AU* + S )dx, 


( 7 ) 


{ ~a(u i+ 1 - Ui) + u(p i+ 1 - p^ 
Y [ u i + 1 ~ «*)] 



( 8 ) 


The first term of the Eq. (8) is the result of the exact integration of -AUj, over the cell C. The integration 
of the second term is not exact and therefore is the source of the overall discretization error. We will further 
discuss this discretization error in the following sections. 

We construct the spatial discretization of Eq. (4) by distributing the cell residuals to the nodes: 


^ = -f(S+$ L + £"$*), (9) 

clt hi 

where <& L and denote the cell-residuals over the left and right cells of the node i, respectively, and hi is 
the dual volume (see Fig. 1) defined by 


hi = 


hL + hR 
2 


hL Xi %i—li hR — Xi. 


(10) 


Note that the pseudo time derivative on the left hand side is retained here just for the sake of illustration. It 
will be ignored in the implicit formulation that follows. Our aim is to directly solve it for the pseudo steady 
state, which corresponds to the next physical time. In this sense, our method is not a dual-time stepping 
method. The distribution matrices B~ and B + (Fig. 2) typically do not affect the order of accuracy of 
the discretization scheme; it is determined by the cell-residuals, and therefore our focus will be on the 
residual evaluation. (Refer to Ref. 1 for formulation of the distribution matrices B~ and B + ). 
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5 <D ; fi + < D, 
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SO 


+1 


i+1 


i+2 


Figure 2. Residual distribution to the nodes. 


A. Second-Order Discretization 

We may use a simple trapezoidal rule for the integration of the source term S over the cell C and arrive at 
the following cell residuals: 


<j> c 


-a(ui+ 1 - Ui) + v(pi + 1 - Pi) - (xi+i 




a 

At 


(^£+1 H - ' Ui )/^ 


k + 1 


(^£+1 'U'i (•^i+ 1 *^i)(Pi+l "I” Pi ) /^) /-^r 


(*££+1 *££)(<S£+i H - ^i)/2 

0 


n— l,n 


(ii) 


where k and n are the Newton iteration counter (as described below) and the physical time index, respectively. 
Note that the second term of the Eq. (11), which is computed at the two previous physical time steps, is 
constant during the Newton iteration, and thus will not contribute to the Jacobian. 

The implicit formulation is defined by 


U fc+1 = U fc + AU fc , 


( 12 ) 


where U = (ui,pi,v,2,P2, ■ ■ ■ ,un,Pn ) and k is the iteration counter. The correction AU fc = U fc+1 — U fc is 
determined as the solution to the linear system: 


5Res A U fe 

dU 


— Res'b 


(13) 


where Res fc is the right hand side of Equation (9), which is the unsteady residual vector evaluated by U fc . 
Note that the pseudo time derivative has been ignored. The Jacobian matrix is sparse, having three 2x2 
blocks in each row for all interior nodes and two blocks for boundary nodes. The i-th pair of rows of the 
linear system is given by 

J i -iAU-L 1 + JiAU f + J i+ iAUf +1 = -1(S+$ L + £T$*) fe , (14) 

hi 


where AU^_ 1 


(A«ti,Ap?_i), AU? = (Auf.Ap*), A U? +1 = (AuJ +1 , Ap* +1 ), 

1 d(H+$ L ) 
i_1 “ hi dUi- 1 ’ 

1 (d{B+<S> L ) ( dB~<f> R )\ 

1 ~ hi\ dUi + dUi ) ’ 

1 d(B~^ R ) 

!+1 hi dU i+ i 


(15) 

(16) 
(17) 


We note that the derivative of the distribution matrices, B~ and 1? + , are zero for linear advection-diffusion 
problems and for nonlinear problems with constant a/v values, which, for simplicity, are the only nonlin- 
ear system considered in this paper. However, the proposed schemes are applicable to general nonlinear 
advection-diffusion problems. 


5 of 29 


American Institute of Aeronautics and Astronautics 


The residual Jacobians needed in the Newton solver are exactly in the same form as the above equations, 
but the derivatives of the cell-residuals now include the contribution from the physical time derivative: 


d®R 


d$ L 

dU t 


a — Iir—— 
2 At 

—v 

1 

hR 

_ TV 

_ 2T ; 

, a 

—a — n L — — 

2 At 

V 

1 

h L 

TV 

2 T r 


(18) 


(19) 


At each physical time level n, we solve the pseudo steady problem by Newton’s method with the current 
solution at n as the initial solution. This results in second-order discretization, which was demonstrated 
with several examples (linear and nonlinear) in Ref/ 

We now focus on the truncation error of the source term discretization as it is needed for our discussion 
in the next sections where we show two different approaches to obtain higher-order RD schemes. We 
demonstrate the truncation error of the system by considering the residual of the second equation in our 
hyperbolic advection-diffusion system: 


= (u i+ 1 - Ui - h{jp i+ 1 +pi)/2)/T r , 


( 20 ) 


where h is the width of the cell C (i.e. Xj+i — xi). We now expand the around the node i to obtain the 
truncation error ( T.E .) for the second equation, 9 r p, which after some algebra and rearrangements becomes: 


T.E. (d T p) = 

= o(h 2 ) 


h h 2 

( d x Ui Pi) T {d xx Ui d x Pi) T (A 

2 o 


xxx'Ui ^&xxPi) d - 0(h ) 


/T r 


(21) 


where the first two terms will simply vanish as they satisfy the original equation in the pseudo steady state. 
The non-vanishing last term makes the current discretization second-order. The same results are obtained 
with the first equation and therefore is not repeated here. The above analysis shows that the error arises 
from the source term discretization, not from the flux balance term, which is exact in one dimension. It 
implies that we can achieve high-order accuracy only by improving the source term discretization. In fact, 
in Ref., 1 : a fourth-order finite-difference scheme is constructed in this manner, i.e., by a high-order source 
term discretization with explicit pseudo-time stepping targeted for steady problems. 

Note that the above truncation error analysis is based on the cell-residual expanded at a node. Hence, 
the nodal residual has the same order as the cell-residuals because it is constructed as a weighted sum of the 
cell- residuals as shown in Eq. (9). Note also that the truncation error analysis based on the cell-residual is 
local to the cell and thus valid for any size of the cell, i.e., valid for uniform and nonuniform grids. In the 
rest of the paper, the truncation error analysis is all based on the cell-residual. 

In the next sections, we show two different techniques that will lead to higher-order schemes. In partic- 
ular, we discuss on 1) reformulation of the source term with its divergence form, and 2) new source term 
discretization technique which acts as a correction to the trapezoidal rule. 


IV. High-Order RD Scheme with Source Term Reformulation (Third- &; 

Fourth-Order) 

Consider the source term, S , as the divergence of a function f s : = S. We now replace the source 

terms with their divergence forms, and rewrite our first-order hyperbolic advection-diffusion equation as 

d T u = -ad x u + vd x p + d x fS, (22) 

d T p = (d x u — d x fp )/T r , (23) 

where and are the divergence formulation of the source terms in the u and the p equations, respectively. 
With the above reformulation of the advection-diffusion equation, which will be discussed in more details, 


6 of 29 


American Institute of Aeronautics and Astronautics 


the residual evaluations of the system becomes exact with no special discretization scheme: 





(-AV x + f s x )dx, 


(24) 


! -a{u i+ 1 - Ui) + v(p i+ i - Pi) + (fu)i + 1 - {fu)i 

1 • (25) 

— [u i+ i - Ui - (fp)i+ 1 + {fp)i)] 

We remark that even though the residual evaluation of our reformulated advection-diffusion system is exact, 
the overall accuracy of the scheme depends on the accuracy of the divergence formulation of the source terms 
and how it is formulated. 


A. Third-Order Scheme with Divergence Form 

Following the divergence formulation introduced in, we write the source flux in a more general form: 


f S 


m — 3 ( _i 'in— 1 

V 1 \ (x — Xi) n d x n-lS, 

Z ' 71 . 1 


n —1 


(26) 


= (x ~ Xi)S - Xi) 2 d x S +^(x - Xi) 3 d xx S - - Xi) 4 d xxx S + . .., (27) 

where the source term S and its derivatives S x , S xx , etc. are not constants but functions that vary in 
space. We remark that the source flux, f s , is not necessarily a polynomial of order m because it depends 
on the derivatives of the source term S. In this paper, we do not consider high-order schemes that require 
evaluation of third or higher derivatives. The divergence of the source flux, d x f s , for the third-order accuracy 
(i.e. m = 3) is 

d x f s = S+^-d xxx S = S + 0(h 3 ), (28) 

o 


which is identical to the original equation up to the third-order. We remark that this error does not 
necessarily limit the maximum order attained by numerical schemes. In the next section, we will show that 
a fourth-order scheme can be constructed by a slightly modified divergence form with m = 3. 

Discretization of the reformulated hyperbolic system leads to the cell residual, for example, for the second 
equation as 


= [ui+i-Ui-{fp)i+i + (fp)i]/T r , 


(29) 


where special care must be taken when we evaluate ff and ff, :i (because of the presence of aq in the 
divergence formulation) as they depend on whether the cell residual is being distributed to the left or to the 
right node of the cell (Fig. 3): 




[«i+l “ U i - Up)i+ 1 + (fp)i\ / T r, 


h 2 h 3 

^ 2+1 aq hftSi - )_i + 3 X ~^~^ x x^i+i 


/T r , 


[Ui + 1 -Ui - (fp) z + 1 + (fp )i] /T r , 
h? It ^ 

'Ui ~~p~d x Si + ~z~^xx^i 

2 D 


/T r . 


(30) 


(31) 


The source term discretization is, therefore, not conservative, which is natural and appropriate because the 
source term does not have a conservative property and should not be discretized in such a way. If it is 
conservative, the global sum of the cell-residual for the source term will depend only on the boundary data 
(telescoping property), which is wrong for the source term. 

We now show that the truncation error ( T.E .) of the resulting RD scheme with the above divergence 
formulation of the source term (with to = 3) is in fact third-order. Again, for demonstration purposes we 
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o 


i-1 


O c,_1 <J> C ' 




-o 

i+1 


Figure 3. Cell residual evaluation with divergence form of the source term. 


consider the second equation (i.e. d T p)\ the same process can be repeated for the first equation. We first 
expand the source fluxes around the node i: 


^ d x f s dx = ff +1 - ff = h R d x f s + + 0{h 4 ), 


where the d x f s is given by Eq. (28), and 


(32) 

(33) 

(34) 

Using Eq. (29), we evaluate the truncation error of the equation d T p by substituting Eqs. (33) and (34) into 
Eq. (32) and expanding all the terms in Eq. (29) around the node i: 


d ~f s = «.s + fa~s+fw; + o(/.‘), 

S 2 h 3 4 

8 xxx f = 8 XX S 4“ h R d xxx S 4- h R d xxxx S 4- —^~d xxxxx S 4- 0{h ). 


T.E. (d r p) = 


(d x Ui Pi) 4- [pxx'U’i 8 :r p, ) 4~ (d xxx Ui d xx pi) 

2 o 


/T r 


h 3 

~prr(8 xxxx Ui 14 d xxx Pi) 4- 0(h ) 


24 

= 0 4 -0(h 3 ) 


/T r 


1351 


where the first three terms vanish as they satisfy our original equation. The presence of the last term confirms 
that the proposed divergence formula with m = 3 makes our scheme third-order accurate. The same result 
is obtained for the first equation with the divergence formulation of the corresponding source terms, and 
therefore the process is not repeated here. 

We remark that the cost of this new third-order accurate RD hyperbolic advection-diffusion scheme using 
the divergence formulation of the source terms is only the evaluation of the first and second derivatives of 
the source terms; i.e. d x S and d xx S. Note that the source flux is third-order accurate (see Eq. 28). Thus, 
according to the Eqs. (30) and (31), we need second-order and first-order accurate discretization for d x S and 
d xx S, respectively. We derive these discretization by constructing a quadratic function between the nodes 
i and its two neighbors i — 1 and i 4- 1 to arrive (after some algebra) with the following formulas that are 
applicable to the internal nodes of general arbitrary grids (uniform and nonuniform): 


d x Si = b 2 ^ 1 - 1 , + h L] S \ + h L Sl+1 + Q(h 2 ) 


9 XX Sj, — 


h R h L {h R 4- h L ) 

h R Sj-i — (h R 4- hi,)Si 4- hL,Sj + \ 
h R h L (h R + h L )/2 


■0(h), 


(36) 

(37) 


where h R and h R are defined in Eq. (10). The corresponding one-sided formulas for the boundary nodes of 
general arbitrary grids are 


d x S 1 

8 x Sn 

8 X xS i 

8 xx Sn 


-h L .(h L . + 2 h R )S 1 + (h R + h L *) 2 S 2 - h 2 R S 3 


h R h L *(h R 4- Hl*) 
h R (hL + 2h R *)SN — (h R * 4- HlYSn-i 4- h 2 R ,SN-2 
h R -h L (h R • 4- h L ) 

hL*S\ — (h R + hL*)S 2 + h R S 3 n(u\ 
h R h L .(h R + h L .)/2 + 1 

—h R *S]y 4- (h R * + h^jS^-i — h R SM-2 
h R *hL(h R * 4- Hl )/ 2 


0(h 2 ), 


o(h 2 ), 


0(h), 


(38) 

(39) 

(40) 

(41) 
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where hi* and Hr* are defined as 


hi,* =x 3 - x 2 , h R * = x N -i - x N - 2 - (42) 

It is clear from the above formulas that each derivative can be computed in a three-point stencil. Con- 
sequently, the residual at a node is defined in a five-point stencil in the interior, a four-point stencil at the 
nodes adjacent to the boundary, and a three-point stencil at the boundary nodes. In the next section, we 
demonstrate that a fourth-order scheme can be constructed without extending the stencil. 


B. Fourth-Order Scheme with Divergence Form 

Here, we show how a simple modification to the presented divergence form upgrades the scheme order by 
one order; that is our third-order scheme becomes fourth-order with no additional cost. To gain an order, 
we propose the divergence form of the source terms (i.e. function f s ), to be defined as 


m> 3 


/ s = E 


1 (x-x) n d x ~-iS, 


= (x - x)S - ^(x - xfd x S A ^(x - xfd xx S - -^-(x - x) 4 d xxx S + ..., 


(43) 

(44) 


where x = (x^ A Xj+i)/2. Note that the previous divergence form was formed around the x, : while this 
new formulation defines the divergence function around the mid-point x. Similar to the divergence form 
presented in the previous section, the source term and its derivatives are not constants but functions that 
vary in space. And therefore, when the source flux is differentiated it recovers the original source term up 
to 0(h m ) around the node i. For m = 3, we have: 


d x f s = S+ {x R X)3 d xxx S = S + 0(h 3 ). (45) 

o 

Again, it does not limit the maximum order of accuracy for numerical schemes, and it is possible to construct 
even higher-order schemes. 

We now show that this third-order source divergence formulation will result in a fourth-order accurate 
RD scheme. We do this similarly to the process explained in the previous section except that the second 
and third derivatives of the f s are now defined as: 


d xx f s = d x s+ {x 2 x)2 d xxx s+ {x 6 x)3 d xxxx s + o(h 4 ), 


&xxxf — & XX S A (x x)d xxx S A (x x) d xxxx S A 


(x — x) c 

6 


d xxxxx S A 0(h 4 ). 


(46) 

(47) 


Again, for discussion purposes, we consider the second equation of the hyperbolic advection-diffusion system, 
d T p. The truncation error of the p equation after substituting Eqs. (46) and (47) into Eq. (32) becomes: 


T.E. (d T p) = 


h h? 

{d x Ui Pi ) A —R-(d X x'Ui ^xPi) A ~^~{O xxx Ui 9 xx Pi ) 
2 6 


/Tr 


h 3 h 4 5 

-^^(^xxxx'U’i & xxxPi ) 4” i2Q iPxxxxx^i xxxxPi ) 4” 0(h ) 


/T r 


= 0 A 0(h 4 ), 


(48) 


which is fourth-order accurate because of the cancellation of the first four terms of the truncation error 
equation. Another great property of this new divergence formulation is the equal evaluation of the ff +1 — f f 
regardless of whether the source flux is being transferred to the node i or i A 1. This greatly simplifies the 
implementation of this form of divergence formulation. 

Note that the identical residual is distributed to the left and right nodes in this scheme. It then appears 
that the source term discretization is conservative. However, it is actually not conservative for the source 
term because the cell-residual for the source term depends on the coordinate of the mid-point of the cell 
and thus it is not telescoping. Again, this is natural and appropriate for the source term, which has no 
conservation property and should be local. 
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V. Higher-Order RD Scheme with Generalized Trapezoidal Rule (Third-, 

Fourth- &c Sixth-order) 


In the previous section, we reformulated the original hyperbolic advection-diffusion system with a gener- 
alized divergence form of the source terms and arrived at third- and fourth-order RD schemes. We specifically 
showed that both the third- and the fourth-order RD schemes developed with the divergence formulation of 
the source terms require the evaluation of the first and second derivatives of the source terms. Fifth- and 
higher-order schemes can be systematically constructed with the divergence formulation, but may require 
the evaluation of third- and higher-order derivatives that would extend the stencil to the neighbors of the 
neighbors and beyond. From a practical point of view, such high-order schemes are not very attractive, and 
thus not considered here. 

In this section, we introduce a different technique to develop even higher-order schemes without the 
need to evaluate gradients beyond the second-derivatives. We also show that with this new technique the 
fourth-order RD scheme is even less computationally intensive than the third- and fourth-order RD schemes 
that are developed with the divergence formulation of the source terms. 

Consider the vector form of our first-order hyperbolic advection-diffusion system: 


<9U 

Ih 


+ A 


<9U 

dx 


= S. 


(49) 


We showed in Section III that the source term discretization with the trapezoidal rule provides a second-order 
accurate scheme. This trapezoidal rule can be written as 



Sdx = -f(S L + Sr), 


(50) 


where for the second-order scheme Sl = Si and Sr = ; i.e. the arithmetic averaging of the source terms 

between the left and the right nodes (Fig. 4). 



Figure 4. Source term discretization. 

Generalizing the trapezoidal rule, we propose the following formula for the left and right source terms, 
Sl and Sr respectively: 

S L - S i + Ct'd x S i + C£d xx S i , (51) 

S R = S i+1 +C?d x S i+1 +C?d xx S i+1 , (52) 

where C\ J and are constants of 0(h), and Cf and C 2 are of 0(h 2 ). A somewhat similar approach is 
taken in Ref., introduced for upgrading the order of accuracy of a finite- volume scheme, where not only the 
source term but also the interface flux need to be modified for high-order accuracy. We find these constants 
by making sure that the nodal residuals of the two equations in our hyperbolic system are accurate up to the 
desired order of accuracy of the scheme we are seeking to develop. With the above equation, the maximum 
possible order of accuracy is six. We also note that the fifth-order RD scheme with the above proposed 
source term discretization becomes mathematically sixth-order. 

A. Third-Order RD Scheme with Generalized Trapezoidal Rule 

A third-order RD scheme can be obtained if the coefficients of the proposed discretization formulation satisfy 
the following equations: 

Cf + Cf = 0, C^h R + C^ + C^ = -^, C?h R + (53) 

which can be derived from the truncation error analysis. Here, there are four unknowns for three equations; 
thus there are infinite combinations with the C ^ as a free parameter. We also note that there is no relation 
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between this third-order RD scheme and the third-order scheme developed with the divergence formulation; 
the nodal residual of the source terms obtained with the previous third-order scheme depends only on the 
information from its immediate node neighbor and not the node itself. On the other hand, the new third- 
order RD scheme that is developed with a correction to the trapezoidal rule depends both on its immediate 
node neighbor and the node itself. Consequently, it is not possible to reproduce the previous third-order 
scheme by any choice of the coefficients. Note, however, that this scheme has the same stencil with the 
previous third-order scheme: a five-point stencil in the interior, a four-point stencil at the nodes adjacent 
to the boundary, and a three-point stencil at the boundary nodes. The same is true for the fourth-order 
scheme derived in the next section. 

Assuming C ^ = —Hr/ 6, the left and right source functions of the source discretization may be chosen 
as: 

S L = Si+bfdxSi^+^d^Si (54) 

S R = S l+1 - ^d x S i+1 - ^-d xx S l+1 . (55) 

o 1000 

We can now evaluate the truncation error of the hyperbolic advection-diffusion system. Here we show the 
truncation error of the second equation; the same will be true for the first equation as well. The cell residual 
of the second equation with the third-order RD scheme with the generalized trapezoidal rule is: 

= ~a(u i+1 - Ui) + u(p i+1 - pi) + ~y(Sl + Sr), 

= -a(u i+ 1 - ui) + v{p i+ 1 - Pi) + + Si) 

- {d x S i+1 -d x Si)--^{d xx S i+1 -d xx Si). (56) 

We then expand the cell residual around the node i to arrive with the following truncation error: 

T.E.(d r Ui) — ( ad x Ui H- vd x pi H- S H- ^ ( ^^xx^i vd X xPi & x Si) 

h 2 

H - — ( Q J &xxx'U J i “I” V^xxxPi “1“ ^xx^i) 

6 

h 3 

+ xxxxPi H - 0.988 &xxx ) + 0(h 4 ) 

= 0 + 0(h 3 ), (57) 

where the first three terms vanish but the last term makes the scheme third-order accurate. As we will show 
next, smaller values for the coefficients C-f or move the schemes to fourth-order accurate. This is further 
shown in the following section. 

B. Fourth- Order RD Scheme with Generalized Trapezoidal Rule 

Fourth-order accuracy is achievable by a simple adjustment to the condition given in Eq. (53): 

C? + C? = 0, C*h R + CL + C* = - h f , °th R + C* = -^ |. (58) 

6 2 12 

Again, the coefficients cannot be determined uniquely, and there are many solutions. A particularly inter- 
esting solution is the following: 

cf = + ^, cf = 0, Cf = -^, Cf = 0. (59) 

As will be shown shortly, it leads to a fourth-order RD scheme with the source term discretization by Eqs. (51) 
and (52): 

S L 

Sr 
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6 

S l+ 1 - ^d x S i+1 , 

n 


( 60 ) 

( 61 ) 
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It is remarkable that fourth-order accuracy is achieved without the second derivatives, which is less expensive 
than the third-order schemes developed in the previous sections. It is, of course, possible to develop a fourth- 
order RD scheme with addition of the second derivatives in the source term discretization. For example the 
following constants will also result in a fourth-order RD scheme: 


C 


L 

1 



ci 


24 


C 


R 


h R r <R _ , h R 
”T’ ° 2 ~ + 24’ 


(62) 


which makes this fourth-order RD scheme identical to the fourth-order RD scheme constructed by the 
divergence formulation in Section B. Although it is possible to construct, these fourth-order schemes are not 
very attractive as they require second derivatives of the source term, and therefore not considered in this 
paper. 

We now prove the fourth-order accuracy of the scheme by evaluating the cell residual of, for example, 
the first equation; i.e. : 


= ~a(u i+ i - Ui) + u(p i+1 - pi) + ~^{S L + Sr), 

h h? 

= -a(u i+1 - Ui) + v(p i+1 -p^ + -^(S i+1 + Si) - ~^(d x S i+ i - d x Si). (63) 

Expanding the cell residual around the node i. we obtain the truncation error of the first equation (after 
some algebra) as 


T .E .(d T Ui) — ( dd x Ui -\~ lsd x Pi “t~ *5)) -t- ^ ( d^xx^i T R^xxPi d - ^x^i) 
h 2 

d" ~ TT ( 0‘dxxx'U'i d” R^xxxPi d” &xx^i) 

6 

h 3 

+ Q’dxxxx'Ri d” vdxxxxPi d” ^xxx^i) 

/l 4 5 

d“ ( R^xxxxx'Ri d" R^xxxxxPi d* R^xxxx^i') d” 0{}P ) 


120 

0 + O(h 4 ), 


6 


(64) 


where the first four terms of the above equation vanish (because of consistency relations). Similar conclusion 
is obtained for the second equation and therefore it is not repeated here. Also, note that the accuracy 
is achieved for general arbitrary grids. The additional cost of upgrading the second-order scheme to the 
fourth-order RD scheme is only due to the evaluation of the second-order accurate first derivative of the 
source terms. The general second-order accurate first derivative formulation for arbitrary grids is provided 
in Sec. A. 


C. Sixth-Order RD Scheme with Generalized Trapezoidal Rule 

In this section, we present a new sixth-order RD scheme with relatively similar cost to our newly introduced 
third-order RD schemes. Specifically, we show that there exist a unique fifth-order RD scheme with the 
generalized trapezoidal rule, and that mathematically becomes sixth-order RD scheme. Seeking fifth-order 
accuracy from the fourth-order RD scheme given in Eq. (58), we find an additional constraint: 


C R h 2 

+ = (65) 

Then, the number of constraints (four) matches the number of unknown coefficients (four). In this case, 
therefore, there exists a unique solution: 


Gf = + 


l R 


C L - —I — h R 
° 2 " + 60’ 


Cf = — 


C? = 


60 ' 


( 66 ) 


Interestingly, the above unique coefficients satisfy the following constraint as well, which is the requirement 
for sixth-order accuracy: 


cl 

4 


h R + C? 


h 2 r 
30 ' 


( 67 ) 
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Thus, we have developed a sixth-order RD scheme with the generalized trapezoidal rule. We remark that 
these coefficients, unlike the ones proposed for the third-order and fourth-order schemes, are unique. With 
the proposed constants for the sixth-order RD scheme, the source term evaluations at the nodes i and i + 1 
become: 


Sl = Si + ^dzSi + ^d^Si, (68) 

5 oU 

S R = S i+1 - f ^fd x S i+1 + f ^d xx S i+1 . (69) 

Note that it requires only the first and second derivatives of the source term, and no higher-order derivatives 
are required. 

Similar to the process explained for the fourth-order scheme, we prove the order of accuracy of the scheme 
using the above proposed source discretization by first evaluating the cell residual of the hyperbolic system. 
Consider the cell residual of, for example, the first equation ($„ ): 

= -a(ui + i-Ui) + u(p i+1 -p i ) + ^Y(S L + S R ), 

= -a(u i+ 1 - Ui ) + v{p i+ 1 - p^ 

+ + Si) — ^(t? x Si+i — d x Si) + Y^(dxocSi + 1 + d xx Si). (70) 

Expanding the cell residual around the node i. we obtain the truncation error of the first equation (after 
some algebra) as 


T.E.(d T Ui) = 

hn 

(- ad x Ui + vd x pi + Si) + — ( -ad xx Ui + vd xx pi + d x Si) 

+ 

h 2 

— — ( Q'dxxx'U'i vdxxxPi H - &xx&i) 
6 

+ 

h 3 

( Gjdxxxx'U'i H - vdxxxxPi H - ^xxx^i) 

+ 

-j^2q ( Q'^xxxxxW'i H - vdxxxxxPi H” &xxx x&i) 

+ 

h R 

Y2q ( Q'^xxxxxx^i “1“ vdxxxxxxPi ^xxxxx^i) 

/ijj 21 

T r 040 ^ Q’&xxxxxxx'Ui T ^dxxxxxxxPi T ^xxxxxx-Si) T 0(h ) 

= 0 + 0(h 6 ). (71) 

The similar result is obtained for the second equation, and therefore is not repeated here. The above 
truncation error analysis reveals that this powerful sixth-order RD scheme could, in practice, produce almost 


seventh-order accurate results, if the sixth derivative of the source term becomes small; the sixth-order term is 
only due to the presence of (h^/100800)d xxxxxx Si in the last term of Eq. (71). We will show such interesting 
results in the following section. 

We emphasize that the sixth-order RD scheme, similar to the third-order RD scheme, requires the 
evaluation of the first and second derivatives of the source terms. However, these derivatives are now 
required to be evaluated with higher-order accuracy. For this sixth-order RD scheme, we need third-order 
and second-order accurate evaluations of the first and the second derivatives of the source terms, respectively, 
which can be obtained by a cubic fit. This makes the developed sixth-order scheme slightly more expensive 
than the third-order scheme. Nevertheless, the proposed sixth-order scheme is exceptionally simple and 
affordable. 


VI. Results 

In this section we present the results in three categories: 1) steady advection-diffusion equation for high 
Reynolds (or Peclet) number applications, 2) unsteady linear advection-diffusion, and 3) unsteady nonlinear 
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advection-diffusion problems. We obtain the results with all the proposed RD schemes and compare them 
with the second-order RD scheme. The order of accuracy results are also compared and presented for each 
example. 

For all cases, we employ the Newton-GS solver as described in Section III. It is essentially Newton’s 
method for the second-order scheme, but an approximate Newton method for higher-order schemes because 
the Jacobian matrix derived from the second-order scheme is used for all higher-order schemes. For unsteady 
problems, the same solver is used to solve the implicit-residual equations or equivalently to compute the 
pseudo steady solution at every physical time step. 

A. Steady Linear Advection-Diffusion 

Consider the advection-diffusion equation in a; £ (0, 1) with u( 0) = 0 and m(1) = 0: 

dtu + ad x u = vd xx u + s(x), (72) 

where 

s(x) = v-K 2 s\h(ttx) + a7r cos(7ra). (73) 

The above problem has a fixed analytical solution of sin(7ra;) for any advection and diffusion coefficients. 
We solved this problem using the hyperbolic advection-diffusion method discussed in Section II with the 
proposed high-order RD schemes (see Fig. 5. Note that we intentionally used a very coarse grid to show 
that the high-order schemes could produce a much more accurate solution on such a “bad” grid.) We 



(a) u{x) 



(b) p(x) = du/dx 


Figure 5. Comparisons between the exact and the numerical results using the proposed high-order RD schemes 
for the steady linear advection-diffusion problem {Re = 1) on a nonuniform grid with N = 10. 

used ranges of nonuniform grids and solved the problem with the Newton-GS method. For each Newton 
iteration, the GS relaxation were conducted until three orders of magnitudes reduction is achieved for the 
linear system. The computations were continued until the residuals of both the solution u and the solution 
gradient p were reduced by eight orders of magnitude. The convergence results are shown in Table 1, where 
the convergence of the GS relaxation is clearly 0(N ) and is independent of the order of accuracy of the RD 
scheme. Also, the solution was converged with the same number of GS relaxations and Newton iterations 
regardless of the RD scheme order. Note that the solutions were converged with a very small number of 
Newton iterations, typically less than ten; this is exceptionally remarkable for the approximate Jacobian 
(second-order) formulation for higher-order schemes. 

The accuracy of the proposed RD schemes were verified by computing the Li = Y^h=i(Ui Xact ~ Ui)/N. 
These results are shown in Fig. 6 for the third-, fourth-, and sixth-order hyperbolic RD schemes (see also 
Tables 2, 3, and 4). The RD schemes developed with the divergence formulation are denoted by RD-D , while 
the schemes developed with the generalized trapezoidal rule are denoted with RD-GT. The results show that 
all the RD schemes achieve the design order of accuracy for both the solution u and the solution gradient 
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Table 1. Steady linear advection-diffusion problem, Re = 1 (Convergence criteria: Residuals <10 8 .) 


Number of nodes 

RD Scheme Order 

GS relaxations/Newton iteration 

Newton iteration 



High-Order Technique 

High- Order Technique 



RD-D 

RD-GT 

RD-D 

RD-GT 


3rd 

168 

169 

10 

10 

50 

4th 

168 

168 

10 

10 


6th 

- 

168 

- 

10 


3rd 

325 

325 

8 

8 

100 

4th 

325 

325 

8 

8 


6th 

- 

325 

- 

8 


3rd 

670 

670 

7 

7 

200 

4th 

670 

670 

7 

7 


6th 

- 

670 

- 

7 


3rd 

1015 

1015 

7 

7 

300 

4th 

1015 

1015 

7 

7 


6th 

- 

1015 

- 

7 


3rd 

1703 

1703 

7 

7 

500 

4th 

1703 

1703 

7 

7 


6th 

- 

1703 

- 

7 


6rd 

3416 

3416 

7 

7 

1000 

4th 

3416 

3416 

7 

7 


6th 

- 

3416 

- 

7 


Third-Order RD Schemes 


Fourth-Order RD Schemes 


Sixth-Order RD Scheme 

H-U (RD-GT) 
q p=U A (RD-GT) 


H-u (RD-GT) 
A p=u x (RD-Gf) 


H-U (RD-GT) 
_ 0 _p=u x (RD-Gf) . 

4 

4 

A u (RD-D) 


A u (RD-D) 



^ p=u x (RD-D) 


^ p=u x (RD-D) 










, - 

J- 

'«(' 



f 

Sl„p C 2..." 

go 2 ‘ 4 

slop ' 2 — 


-6 

slop ' 3 

-6 

Slope 3 . - • 7/ 


-* 

Slope * 

-« 

Slope 4 


-10 


-10 






Slope 6 Aj-' 







5 -3 


-3 

-3 -2.5 -2 -1.5 -1 -0 

log io h 


(a) third-order 


(b) fourth-order 


(c) sixth-order 


Figure 6. L\ error of the proposed high-order RD hyperbolic schemes for the steady linear advection-diffusion 
problem. 
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p. It can be seen also that the difference between the two versions of third-order RD schemes is minor, and 
the same is true for the fourth-order RD schemes. 


Table 2. Spatial accuracy with the third-order RD-GT scheme for the steady linear advect ion-diffusion problem 

(Re = 1). 


Number of nodes 

Li error of u 

Order 

L i error of p 

Order 

10 

9.54E-03 


8.59E-02 


25 

4.63E-04 

3.30 

4.41E-03 

3.24 

50 

4.85E-05 

3.25 

4.63E-04 

3.25 

100 

5.46E-06 

3.15 

5.16E-05 

3.17 

200 

6.46E-07 

3.08 

6.04E-06 

3.09 

300 

1.88E-07 

3.04 

1.75E-06 

3.05 

500 

3.99E-08 

3.03 

3.71E-07 

3.04 

1000 

4.93E-09 

3.02 

4.57E-08 

3.02 


Table 3. Spatial accuracy with the fourth-order RD-D scheme for the steady linear advection-diffusion problem, 
Re = 1. (Note: fourth-order RD-D and RD-GT schemes produce almost identical result.) 


Number of nodes 

L\ error of u 

Order 

L i error of p 

Order 

10 

8.38E-03 


6.36E-02 


25 

2.84E-04 

3.69 

1.82E-03 

3.88 

50 

1.74E-05 

4.03 

1.02E-04 

4.16 

100 

1.04E-06 

4.07 

5.77E-06 

4.15 

200 

6.23E-08 

4.06 

3.35E-07 

4.11 

300 

1.21E-08 

4.04 

6.44E-08 

4.07 

500 

1.55E-09 

4.02 

8.15E-09 

4.05 

1000 

9.30E-11 

4.06 

4.94E-10 

4.04 


B. Steady Boundary Layer Problem 

Consider the advection-diffusion equation in x £ (0, 1) with u( 0) = 0 and m(1) = 1: 

d t u + ad x u = vd xx u + s(x), (74) 


where 


7T 

s(x) = — [acos(7ra) + nu sin(7r:r)] , Re = a/v. 


Re 


(75) 


This is a boundary layer problem with a non-trivial steady state solution in the diffusion limit as a result of 
the source term addition. This equation develops a very narrow boundary layer near the right boundary 
( x = 1) when the advection term becomes dominant. The exact steady state solution to this problem is 
given by (see also Fig. 7). 

p -Re _ ( x-l)Re i 

u e * act (x) = 0 -Re_ l + 777 sin (™)- ( 76 ) 


Re 


We chose various Re values ranging from 1 to 10 6 and solved the equation on nonuniform grid sizes 
up to 100000 nodes. Like the previous problem, the solutions were obtained with the Newton-GS method. 
Within each Newton iteration, the GS relaxation were conducted until three orders of magnitude reduction is 
achieved for the linear system, and the residuals of both the solution and the solution gradient were reduced 
by eight orders of magnitude. The convergence data are given in Table 5 for the proposed high-order 
RD schemes. The data are shown for the high-order schemes developed with the generalized trapezoidal 
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Table 4. Spatial accuracy with the sixth-order RD-GT scheme for the steady linear advection-diffusion problem 
(Re = 1). 


Number of nodes 

Li error of u 

Order 

L\ error of p 

Order 

10 

2.60E-03 


1.75E-02 


25 

3.34E-05 

4.75 

1.98E-04 

4.89 

50 

4.87E-07 

6.10 

2.87E-06 

6.11 

100 

4.32E-09 

6.82 

2.65E-08 

6.76 

200 

3.09E-11 

7.13 

1.78E-10 

7.22 

300 

5.80E-12 

4.13 

3.39E-11 

4.10 

500 

3.60E-12 

0.93 

1.97E-11 

1.06 


Re=10 


Re=10 




(a) u(x) 


(b) p(x) = du/dx 


Figure 7. Comparisons between the exact and the numerical results using the proposed high-order RD schemes 
for the steady boundary layer problem (Re = 10) on a nonuniform grid with N = 10. 
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rule. Similar results were obtained with the RD schemes developed with the divergence formulation, and 
therefore not shown. But in general, the high-order RD schemes with the generalized trapezoidal rule 
perform better than the ones developed with the divergence formulation. The latter third-order RD-D 
scheme not only produce slightly larger errors as shown in Fig. 6, but also encounter some convergence 
difficulties (particularly with time-dependent problems); the convergence problem seems originated from the 
lack of diagonal dominance caused by the particular structure of the source term discretization as shown 
in Section A. We will therefore mainly focus on the RD-GT schemes for unsteady cases. The 0(N ) linear 


Table 5. Steady boundary layer problem (Convergence criteria: Residuals < 10 8 .) 


logio Re 

Number of nodes 

RD-GT Scheme Order 

GS relaxations/Newton iteration 

Newton iteration 



3rd 

163 

8 

0 

50 

4th 

163 

8 



6th 





3rd 

324 

7 

0 

100 

4th 

324 

7 



6th 





3rd 

1647 

7 

0 

500 

4th 

1647 

7 



6th 





3rd 

178 

7 

1 

100 

4th 

178 

7 



6th 





3rd 

44 

7 

2 

100 

4th 

44 

7 



6th 





3rd 

42 

8 

3 

500 

4th 

43 

7 



6th 





3rd 

18 

9 

4 

1000 

4th 

19 

7 



6th 





3rd 

24 

9 

5 

10000 

4th 

21 

7 



6th 





3rd 

18 

9 

6 

100000 

4th 

19 

7 



6th 




dependency of the GS relaxations on the grid size is also demonstrated, which is a consequence of solving the 
advection-diffusion equation as a hyperbolic system. We emphasize that this is remarkable because the linear 
convergence is retained for any irregular grid in any dimensions (N is approximately the number of nodes 
in each coordinate direction in two and three dimensions) as demonstrated in Refs. “ It leads to orders 
of magnitude faster convergence in comparison with conventional methods whose convergence is typically 
0(N 2 ) as discussed also in the previous paper. Furthermore, the proposed high-order RD schemes are 
extremely efficient as the solutions were obtained with only a small number of Newton iterations: less than 
10 iterations to reduce the residual by eight orders of magnitude. Moreover, the number of GS relaxations 
and Newton iterations are essentially independent of the scheme order. Considering the fact that the cost of 
one GS relaxation is significantly cheaper than one Newton iteration, we find that the developed high-order 
RD schemes are extremely powerful and efficient. Finally, as in the previous work, we remark that the 
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high-Re cases required extremely fine grids to meet the well-known requirement on the mesh Reynolds- 
number. J If desired, the computations can be performed on substantially coarser grids with more aggressive 
and customized grid stretching. However, we simply refined the grid to meet the mesh Reynolds-number 
requirement because our method is powerful enough to solve the problem very efficiently (i.e. , less than 10 
Newton iterations) even on such dense grids for all third-, fourth-, and sixth-order schemes. The ability to 
efficiently solve the problem on highly refined grids is a great advantage of these schemes. 

The order of accuracy of the proposed RD schemes were also verified for this problem. Figures 8 shows 
the L\ error convergence results, where h is the representative mesh spacing defined by h = /(N — 1). For 
discussion purposes, we present the accuracy plots for Re = 1 and Re = 10; similar results were obtained for 
other Reynolds numbers. These results verify the order of accuracy of the proposed high-order RD schemes 
(i.e. third-, fourth-, and sixth-order) for all the variables and the gradients at all the grid nodes including 
the boundary nodes (see also Tables 6, 7, and 8). Note that when the differences between the numerical 
results and exact values approach the machine zero, the L\ slope approaches zeroth-order. 
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Figure 8. L\ error of the proposed high-order RD hyperbolic schemes for the boundary layer problem on 
nonuniform grids. 


Table 6. Spatial accuracy with the third-order RD-GT scheme for the boundary layer problem with Re = 1.0. 


Number of nodes 

Li error of u 

Order 

Li error of p 

Order 

10 

9.65E-03 


8.68E-02 


25 

4.70E-04 

3.30 

4.47E-03 

3.24 

50 

4.94E-05 

3.25 

4.70E-04 

3.25 

100 

5.56E-06 

3.15 

5.26E-05 

3.16 

200 

6.58E-07 

3.08 

6.16E-06 

3.09 

300 

1.91E-07 

3.06 

1.79E-06 

3.05 

500 

4.07E-08 

3.03 

3.79E-07 

3.04 

1000 

5.04E-09 

3.01 

4.66E-08 

3.02 

2000 

6.34E-10 

2.99 

5.79E-09 

3.01 

5000 

4.89E-11 

2.80 

3.68E-10 

3.00 


C. Unsteady Linear Advection-Diffusion 

Consider the time-dependent advection-diffusion equation in a; £ (0, 1) 


d t u + ad x u = vd xx u. (77) 

The above equation with the following initial condition: 

u(x,t = 0) = sin(Kx), (78) 
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Table 7. Spatial accuracy with the fourth-order RD-GT schemes for the boundary layer problem with Re = 1.0 
Note that when the differences between the numerical results and exact values approach the machine zero 
the L\ slope approaches zeroth-order. 


Number of nodes 

Li error of u 

Order 

L i error of p 

Order 

10 

8.48E-03 


6.43E-02 


25 

2.88E-04 

3.69 

1.84E-03 

3.88 

50 

1.78E-05 

4.02 

1.04E-04 

4.14 

100 

1.06E-06 

4.07 

5.87E-06 

4.15 

200 

6.41E-08 

4.05 

3.42E-07 

4.10 

300 

1.25E-08 

4.04 

6.57E-08 

4.07 

500 

1.61E-09 

4.01 

4.35E-09 

4.05 

1000 

1.06E-10 

3.92 

4.98E-10 

4.06 

2000 

1.42E-11 

2.90 

1.42E-11 

3.40 

5000 

9.69E-12 

0.42 

4.75E-11 

0.01 


Table 8. Spatial accuracy with the sixth-order RD-GT scheme for the boundary layer problem with Re = 10 


Number of nodes 

Li error of u 

Order 

L\ error of p 

Order 

10 

4.81E-02 


5.54E-01 


25 

6.89E-05 

4.75 

6.19E-04 

4.89 

50 

4.54E-07 

6.11 

3.76E-06 

6.11 

100 

3.97E-09 

6.78 

3.21E-08 

6.73 

200 

3.56E-11 

7.59 

4.28E-10 

7.37 

300 

1.14E-11 

3.72 

1.40E-10 

4.23 

500 

9.23E-12 

0.05 

1.19E-10 

0.30 
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where k is an arbitrary constant, has the following exact solution with a periodic boundary condition: 


u exact (x, t ) = e~ K vt sin(n(x — at)). 

A non-periodic solution also exists for the following oscillatory boundary conditions 

u(0,t) = 0, 

u(l,t) = U cos(wf), 


( 79 ) 


(80) 

(81) 


where U is the amplitude of the oscillation and ui is the frequency of the oscillation on the right boundary. 
The exact solution is given by 


u exact {x,t) = Real 




Al,2 — 


a ± \/ a 2 + 4 iuv 
2v 


(82) 


where i = \J — 1. 

We solved the above two problems with the first-order hyperbolic advection-diffusions equation given as 
in Eq. (2) and (3). For each physical time, we reduced the residuals by two orders of magnitude before 
advancing in time. During each time step, we also relaxed the linear system using GS relaxations until two 
orders of magnitude reduction in the linear system residuals was achieved (see Figs. 9 and 10.) We also note 
that more residual reduction in pseudo time may be necessary on more complex problems but two orders of 
magnitude reduction in the residuals were sufficient for the problems presented here. 




(a) u(x,t = 1.0) (b) p(x,t = 1.0) = du/dx 


Figure 9. Time-dependent linear advection-diffusion problem (Re = 33.33) with periodic boundary condition 
on IV = 10 uniform nodes (At = 0.01.) 

We examined the convergence rate of these problems on several uniform and nonuniform grid systems. 
Given in Tables 9 and 10 are the average numbers of GS relaxations per Newton iteration obtained over 
100 time steps for the periodic and the oscillatory boundary condition problems, respectively. Clearly the 
convergence rate of the GS relaxation is of 0(N ), not 0(N 2 ) as typical for numerical methods for the 
advection-diffusion equation. Observe that for most grid systems only two Newton iterations were sufficient 
to obtain accurate solutions regardless of the scheme order. We also note that we used At = 0.01 for all 
the grid systems. The time step is orders-of-magnitude larger than that required for conventional explicit 
schemes, which is limited by 0(h 2 ). Of course, conventional implicit schemes also allow unconditionally 
large time steps, but it requires 0{N 2 ) convergence in an iterative linear solver and potentially a much 
larger number of outer iterations as well if the exact linearization is not possible and Newton’s method 
cannot be constructed. Hence, the method developed here has two major advantages over conventional 
methods: the second order Jacobian formulation and 0(N) iterative convergence in the linear solver. The 
latter advantage can be potentially huge with increase of the grid system as the speed-up factor is O(N) and 
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(a) u(x,t = 1.0) (b) p(x,t = 1.0) = du/dx 


Figure 10. Time-dependent linear advection-diffusion problem with oscillatory boundary condition (u = 77r/2, 
ac = 2-7T, U = 1, v = 1) on N = 10 uniform nodes (At = 0.01.) 


Table 9. Unsteady linear advection-diffusion problem with periodic BC (a = l,i/ = 0.03) on uniform grids. 
Average data over 100 time steps are given. (Convergence criteria: Residuals < 10 — 2 ) 


Number of nodes 

RD-GT Scheme Order 

GS relaxations/Newton iteration 

Newton iteration 


3rd 

5 

4 

25 

4th 

5 

4 


6th 

6 

5 


3rd 

24 

2 

100 

4th 

24 

2 


6th 

24 

2 


3rd 

33 

2 

300 

4th 

35 

2 


6th 

35 

2 


3rd 

55 

2 

500 

4th 

55 

2 


6th 

55 

2 


3rd 

116 

2 

1000 

4th 

116 

2 


6th 

116 

2 
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Table 10. Unsteady linear advection-diffusion problem with oscillatory BC (uj = 77r/2,a = 1.) on nonuniform 
grids. Average data over 100 time steps are given. (Convergence criteria: Residuals < 10~ 2 ) 


Number of nodes 

RD-GT Scheme Order 

GS relaxations/Newton iteration 

Newton iteration 


3rd 

35 

4 

25 

4th 

37 

4 


6th 

45 

4 


3rd 

69 

4 

50 

4th 

72 

4 


6th 

72 

4 


3rd 

136 

4 

100 

4th 

142 

4 


6th 

142 

4 


3rd 

650 

4 

500 

4th 

680 

4 


6th 

680 

4 


3rd 

1283 

4 

1000 

4th 

1332 

4 


6th 

1332 

4 


thus grows for finer grids. Note also that the second-order Jacobian formulation is the advantage of the RD 
scheme over finite volume schemes where the compact Jacobian formulation is only first-order. The results 
also show that the convergence rate is the same for all the developed high-order RD schemes. It means that 
the only cost to these higher order schemes is the evaluation of the first and the second derivatives of the 
source terms and in the case of the fourth-order scheme with the generalized trapezoidal rule, the second 
derivatives are not required. 

We verified the order of accuracy of the proposed schemes on time-dependent linear problems with 
consistent space-time discretization; i.e. third-order with BDF3, fourth-order with BDF4, and sixth-order 
with BDF6. The same spatial order of accuracy can be observed with the A-Stable BDF2 with small 
enough time steps such that the temporal error is comparable to the spatial error (i.e. At 2 ~ h m , where 
m is the scheme order) . But the consistent pair of BDF and spatial discretization allows about two orders 
of magnitude reduction in the number of time steps for the finest grid cases. Note that time-accurate 
computations are started by BDF1 in the first time step with extremely small time step (e.g. At = 10 -8 ), 
and then by higher order BDFs thereafter with much larger time steps (e.g. At = 0.01). This will ensure 
the order of accuracy of the developed higher order schemes through all times. We remark that explicit time 
stepping is not available for time- accurate computations with the hyperbolic system method (see Ref. for 
more details.) 

Figure 11 shows the L\ error convergence for the above time-dependent linear problem with the oscillatory 
boundary condition (see also Tables 11, 12, and 13), where clearly the order of accuracy of the proposed 
schemes are verified for the linear time-dependent problems. The results were obtained at t = 1.0. We used 
the same At among the third-, fourth- and sixth-order schemes and were able to obtain the desired order of 
accuracy. (Note that we showed the RD schemes that were developed with the generalized trapezoidal rule 
as these are more efficient and less expensive than the ones developed with divergence formulation of the 
source terms.) 

D. Unsteady Nonlinear Advection-Diffusion 

Consider the unsteady nonlinear viscous Burgers equation with an unsteady time-dependent source term: 

d t u + d x f = d x (vu x ) + S(x,t), a; €(0,1), (83) 
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Third-Order RD Scheme 


Fourth-Order RD Schemes 


Sixth-Order RD Scheme 





(a) third-order + BDF3 (b) fourth-order + BDF4 (c) sixth-order + BDF6 

Figure 11. Spatial accuracy of the proposed high-order RD-GT hyperbolic schemes for the time-dependent 
linear problem with oscillatory BC on uniform grids. 


Table 11. Spatial accuracy for the linear time dependent problem with oscillatory BC (lj = 7ir/2 ,a = 1.) using 
the third-order RD-GT scheme with the BDF3 time discretization. 


Number of nodes 

At (BDF3) 

L\ error of u 

Order 

L i error of p 

Order 

10 

2.50E-03 

1.40E-04 


2.86E-04 


20 

2.50E-03 

9.06E-06 

3.95 

1.90E-05 

3.91 

50 

1.25E-03 

3.54E-07 

3.54 

7.38E-07 

3.54 

100 

5.00E-04 

2.57E-08 

3.78 

4.71E-08 

3.97 

200 

2.50E-04 

2.59E-09 

3.31 

4.40E-09 

3.42 


Table 12. Spatial accuracy for the linear time dependent problem with oscillatory BC (cj = 7tt/ 2 ,a = 1.) using 
the fourth-order RD-GT scheme with the BDF4 time discretization. 


Number of nodes 

At (BDF4) 

L\ error of u 

Order 

L i error of p 

Order 

10 

2.50E-03 

1.32E-04 


2.83E-04 


20 

2.50E-03 

7.22E-06 

4.19 

1.63E-05 

4.12 

50 

1.25E-03 

1.70E-07 

4.09 

4.01E-07 

4.04 

100 

5.00E-04 

1.04E-08 

4.03 

2.51E-08 

4.00 

200 

2.50E-04 

6.78E-10 

3.94 

1.64E-09 

3.94 


Table 13. Spatial accuracy for the linear time dependent problem with oscillatory BC (uj = 77r/2,a = 1.) using 
the sixth-order RD-GT scheme with the BDF6 time discretization. 


Number of nodes 

At (BDF6) 

L\ error of u 

Order 

L i error of p 

Order 

10 

2.50E-03 

3.77E-06 


2.29E-05 


20 

2.50E-03 

7.19E-08 

5.71 

3.27E-07 

6.13 

50 

1.25E-03 

5.81E-09 

6.20 

3.70E-08 

5.37 

100 

5.00E-04 

3.62E-10 

5.43 

1.34E-09 

6.50 


24 of 29 


American Institute of Aeronautics and Astronautics 


where / = it 2 / 2, v = u, and 


(84) 


S(x,t) = d t u e + -d x (u e ) 2 - d x {u e p e ), 

where p e = d x u e . The source term has been generated by the following function: 


u'(x,t) = Real ( sillh(l /g Z; ) ce‘",) +C, C > 1, 

\ sinh(-\/iw/i/) / 

so that it is the exact solution to Eq. (83) with the boundary conditions defined as 


u(0,t) = C , 

u(l,t) = C + (7cos(wf), 


(85) 


( 86 ) 

(87) 


where ui is the frequency of the oscillation on the right boundary, and U is the amplitude of the oscillation. 
We note that the constant C must be greater than 1 in order for the diffusion coefficient to be positive. We 
solved this time-dependent nonlinear advection-diffusion equation with the following equivalent first-order 
hyperbolic system (see Ref. for more details): 


d T u + d x (w 2 ) = d x p — d t u + S(x,t), 

(88) 

T 

— d T p = (d x u-p/v). 

(89) 


For this nonlinear unsteady problem, the manufactured source term contains terms that are already in the 
divergence form; i.e. the residual evaluation of these terms are exact. The dtu e term is the only term in the 
manufactured source term with a non-exact residual evaluation. The BDF discretization of the dtu term 
in Eq. (88) will not be in the divergence form and therefore will not have an exact residual evaluation. In 
addition, the second equation has a nonlinear source term, p/ v (note that here v = u). We obtained the 
high-order results by evaluating all of these non-exact residuals using the proposed techniques. Newton 
iterations are taken to be converged when the overall residuals are dropped by eight orders of magnitude. 
For each Newton iteration, we relaxed the linear system until the residuals are reduced by two orders of 
magnitude. 

The O(N) convergence rate of the GS relaxations was once again achieved for the time-dependent nonlin- 
ear hyperbolic advection-diffusion system. This is given in Table 14, where the average number of iterations 
were obtained over 1000 time steps (over 17 periods). Note also that the high-order RD schemes converged 
with only ten Newton iterations with the compact second-order Jacobian formulation. It is remarkable that 
such a rapid convergence is achieved for all high-order schemes with the second-order Jacobian. 


Table 14. Unsteady nonlinear advection-diffusion problem (u = a = u) with oscillatory BC (U = 1, C = 2, 
lo = 7tt/2) on nonuniform grids. Average data over 1000 time steps are given. (Convergence criteria: GS 
Relaxation < 1CU 2 ; Newton residuals < 10~ s ). 


Number of nodes 

RD-GT Scheme Order 

GS relaxations/Newton iteration 

Newton iteration 


3rd 

435 

10 

50 

4th 

430 

10 


6th 

431 

10 


3rd 

879 

10 

100 

4th 

868 

10 


6th 

864 

10 


3rd 

1772 

10 

200 

4th 

1749 

10 


6th 

1737 

10 


Shown in Fig. 12 are the order of accuracy plots obtained for this unsteady advection-diffusion problem 
on series of nonuniform grids (see also Tables 15, 16, and 17). The results confirm the high-order accuracy 
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of developed RD schemes for unsteady nonlinear problems on nonuniform grids. We also observe that our 
proposed sixth-order RD scheme, as discussed in the previous section, have in fact produced seventh-order 
accurate results demonstrating its power characteristics in efficiently providing very high accurate solutions 
and gradients (see Table 17 for more details). We note, however, that the seventh-order accuracy is not a 
general feature of the sixth-order scheme. Although it is not clear from the exact solution, it could be due 
to the sixth-order derivative of the source term in the leading truncation error being vanishingly small. 


Third-Order RD Scheme 


Fourth-Order RD Schemes 


Sixth-Order RD Scheme 



-e-p(= u x> 



(a) third-order + BDF3 


lo S.o h 

(b) fourth-order + BDF4 



(c) sixth-order + BDF6 


Figure 12. Spatial accuracy of the proposed high-order RD-GT hyperbolic schemes for the time-dependent 
nonlinear problem (y = a = u) with oscillatory BC ( U = 1, C = 2, oj = 77r/2) on nonuniform grids. 


Table 15. Spatial accuracy for the time-dependent nonlinear problem (y — a = u) with oscillatory BC ( U = 1, 
C = 2, uj = 7tt/2) using the third-order RD-GT scheme and BDF3 time discretization on nonuniform grids. 


Number of nodes 

At (BDF3) 

L\ error of u 

Order 

L\ error of p 

Order 

25 

1.25E-03 

6.83E-05 


3.59E-04 


50 

1.25E-03 

3.68E-06 

4.21 

1.79E-05 

4.33 

100 

1.00E-03 

2.33E-07 

3.98 

1.12E-06 

4.00 

200 

5.00E-04 

1.70E-08 

3.78 

8.23E-08 

3.77 

500 

2.50E-04 

9.58E-10 

3.14 

5.03E-09 

3.05 


Table 16. Spatial accuracy for the time-dependent nonlinear problem (y = a = u) with oscillatory BC (U = 1, 
C = 2, to = 77t/ 2) using the fourth-order RD schemes and BDF4 time discretization on nonuniform grids. 


Number of nodes 

At (BDF4) 

L\ error of u 

Order 

L\ error of p 

Order 

25 

1.25E-03 

6.75E-05 


3.49E-04 


50 

1.25E-03 

6.21E-07 

4.27 

2.71E-06 

4.23 

100 

1.00E-03 

1.85E-07 

4.24 

7.68E-07 

4.40 

200 

1.00E-03 

1.07E-08 

4.11 

3.69E-08 

4.38 

500 

5.00E-04 

3.14E-10 

3.85 

6.33E-10 

4.44 
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Table IT. Spatial accuracy for the time-dependent nonlinear problem (v = a = u) with oscillatory BC ( U = 1, 
C = 2, lj = 7tt/2) using the sixth-order RD-GT scheme and BDF6 time discretization on nonuniform grids. 


Number of nodes 

At (BDF6) 

L\ error of u 

Order 

L\ error of p 

Order 

50 

1.25E-03 

7.99E-08 


5.61E-08 


75 

5.00E-04 

5.42E-09 

6.60 

3.69E-09 

6.71 

100 

5.00E-04 

7.54E-10 

6.85 

5.09E-09 

6.91 

200 

5.00E-04 

4.30E-11 

7.06 

2.83E-09 

7.11 


VII. Conclusions 

In this paper, we have developed a series of efficient high-order Residual-Distribution (RD) schemes 
for general advection-diffusion problems. Third- and fourth-order RD schemes were developed with the 
divergence formulation of the source term. These schemes are very economical as they require only the 
evaluation of the first and second derivatives of the source term, and not of the solution. The first and 
second derivatives need to be second and first order accurate, respectively, on arbitrary grids, and they 
are obtained by a compact quadratic fit. Third-, fourth-, and sixth-order RD schemes are also developed 
with a generalized trapezoidal rule. The third-order schemes require the evaluation of the first and second 
derivatives of the source term, which must be second and first order accurate, respectively. On the other 
hand, the fourth-order scheme requires only the second-order accurate gradients, and therefore is the least 
expensive scheme among all the developed high-order schemes in this paper. All third- and fourth-order 
schemes are constructed within a five-point stencil in the interior, a four-point stencil at the nodes adjacent 
to the boundary, and a three-point stencil at the boundary nodes. For the sixth-order scheme, the evaluation 
of the first and second derivatives of the source term is required, and they must be third- and second-order 
accurate, which is achieved by a cubic fit. It results in the stencil extension by one or two nodes in the nodal 
residual. In addition, the analysis of the proposed sixth-order RD scheme as well as the results presented 
here show that this high-order RD scheme can, in some cases, produce seventh-order results on nonuniform 
grids. An implicit steady solver is constructed based on the Jacobian derived from the compact second-order 
RD scheme. It is demonstrated that the solver achieves a rapid convergence like Newton’s method for all 
high-order schemes despite the fact that the Jacobian is not exact. Specifically, it requires only a small 
number of Newton iterations, typically less than 10, for both steady and unsteady problems, even for highly 
refined grids, up to 100000 nodes. We have demonstrated also that all of these high-order schemes are 
capable of producing both high-order accurate solution and gradient on nonuniform grids very efficiently by 
less than ten Newton iterations. 

The study presented in this paper should be of interest to researchers working on finite-volume schemes 
because the RD scheme is known to be equivalent to the finite-volume scheme in one dimension (see for 
example Ref. ). Specifically, a second-order upwind RD scheme is equivalent to the first-order finite- volume 
scheme with a special form of source term discretization. ' 1 It implies that the developed high-order RD 
schemes may also be implemented in the form of first-order finite-volume schemes with special source term 
discretization formulas. The resulting finite-volume schemes will be different from many other finite-volume 
schemes in that they do not require computations of solution gradients. 

The developed high-order schemes could well bring significant improvements to the numerical methods 
for practical problems such as material thermal response calculations of thermal protection systems of at- 
mospheric entry vehicles,' “ and the experimental aeroheating data reduction,' ’ which are based on 
one-dimensional analyses and routinely used in industries (e.g. NASA). A particularly useful scheme would 
be the fourth-order scheme based on the generalized trapezoidal rule (RD-GT) because it requires only 
the second-order accurate gradients of the source term. Application to these practical problems should be 
undertaken and is left as future work. 

Extensions to higher dimensions are highly desired. To extend the developed high-order schemes to higher 
dimensions, it is necessary to develop a high-order quadrature formula for integrating the flux divergence 
term, which is exact in one dimension but not in higher dimensions. For the source term discretization, the 
divergence formulation can be extended relatively straightforwardly while a discretization formula such as the 
generalized trapezoidal rule remains to be found. We also remark that the BDF2 is A-stable but higher-order 
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BDFs are not and may encounter instability in practical problems for complex equations and geometries in 
higher dimensions. Alternative implicit time integration methods, such as space-time methods, may turn 
out to be useful or necessary for such applications. Finally, extensions to more complex nonlinear equations 
such as the Navier-Stokes equations remain as a challenge. For the compressible Navier-Stokes equations, 
the complete eigen-structure of the whole system has yet to be found. The construction of the upwind RD 
scheme based on a single hyperbolic system, as presented in this paper, is a challenge. To overcome the 
difficulty, a simplified approach has been proposed and demonstrated for a finite- volume scheme in Ref., ' 
which is based on the independent treatment of the inviscid and viscous terms. A similar approach may 
become necessary for the RD schemes. 
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